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This paper summarises an investigation of discreteness effects and the continuum limit for time- 
independent, nearly collisionless A-body systems of charged particles interacting via an unscreened 
1/r 2 force, allowing for bulk density distributions corresponding to potentials that admit both 
regular and chaotic orbits. Both for orbit ensembles and for individual orbits, for N — > oo there 
is a smooth convergence towards a continuum limit. At least for moderately large values of N 
discreteness effects are extremely well modeled by Gaussian white noise with energy relaxation time 
iij, and hence diffusion constant D, consistent with the scaling tn oc(A/lnA)i.D, with A the Coulomb 
logarithm and to a natural 'dynamical' time scale, as predicted by a Fokker-Planck description. 
Discreteness effects accelerate emittance growth for an initially localised ensemble of orbits ('clump'). 
However, even allowing for discreteness effects one can distinguish clearly between orbits which, in 
the continuum limit, feel a regular (nonchaotic) potential, so that emittance grows as a power law 
in time, and chaotic orbits, for which the emittance grows exponentially. For sufficiently large 
N, one can implement a clear distinction between two different 'kinds' of chaos acting in A-body 
systems. Short range microchaos, associated with close encounters between individual charges, 
is a generic feature of the A-body problem, giving rise to large positive Lyapunov exponents xn 
which do not decrease with increasing N even if the bulk potential is mtegrable. Alternatively, 
there is the possibility of larger scale macrochaos, characterised typically by somewhat smaller 
Lyapunov exponents xs, which will be present only if the bulk potential admits global stochasticity. 
Conventional computations of the largest Lyapunov exponent provide estimates of xn, leading to 
the oxymoronic conclusion that A-body orbits which look nearly regular and have sharply peaked 
Fourier spectra are 'very chaotic' However, the 'range' of the microchaos is set by the typical 
interparticle spacing which, as N increases, becomes progressively smaller, so that, for sufficiently 
large A, this microchaos, albeit very strong, is largely irrelevant macroscopically. A more careful 
numerical analysis allows one to derive estimates of both xn and xs- 



I. INTRODUCTION AND MOTIVATION 

To what extent can a 'nearly collisionless' A^-body sys- 
tem such as a nonneutral plasma or a charged particle 
beam be modeled by a smooth density distribution and 
a smooth bulk potential, either statistically or at the level 
of individual orbits? In particular, is there a well defined 
N — ► oo continuum limit? The idealisation of a smooth 
potential is extremely convenient, both conceptually and 
computationally. However, as noted, e.g., by beam dy- 
namicists 0, it is not completely clear when - or if - 
such an idealisation is justified. 

Even assuming that there is a well-defined continuum 
limit, how large must N be before the continuum ap- 
proximation is justified? And to what extent can residual 
discreteness effects be modeled in the context of a Fokker- 
Planck description? The conventional Fokker-Planck de- 
scription 2] 3] was formulated to extract statistical prop- 
erties of orbit ensembles and distribution functions over 
long time scales, assuming implicitly that the bulk po- 



tential is regular. To what extent, then, can Langevin 
realisations |^ of a Fokker-Planck equation yield reliable 
information about individual orbits over comparatively 
short time scales, particularly if the orbits are chaotic? 

This is an issue of practical importance for systems like 
high intensity charged-particle beams. In older, low in- 
tensity beams, the contribution of the space charge to the 
total potential is comparatively minor compared to the 
confining magnetic field, but in high intensity beams the 
space charge can become extremely important. To what 
extent, then, is one justified in idealising the space charge 
by a smooth density distribution that generates a smooth 
macroscopic potential? Is it, e.g., really legitimate to ig- 
nore discreteness effects entirely in a bunch comprised of 
~ 10 9 — 10 10 protons with a transverse emittance of a few 
microns? 

The situation is especially suspect in that experience 
with 'nearly collisionless' self-gravitating systems [f| in- 
dicates that the continuum limit must be subtle. One 
believes that, for N — > oo, orbits in an Af-body system 
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will in fact converge towards characteristics in the cor- 
responding smooth potential. However, TV-body orbits 
evolved in a density distribution corresponding to an in- 
tegrable potential are typically 'very chaotic' with large 
positive Lyapunov exponents even though the integrable 
characteristics have vanishing Lyapunov exponent! 

For flows in smooth potentials, it is straightforward to 
distinguish macroscopically between phase mixing asso- 
ciated with regular versus chaotic orbits @|. In particu- 
lar, chaotic flows tend to mix exponentially fast, whereas 
regular flows only mix as a more modest power law in 
time. Will this distinction persist in TV-body systems? 
One might, e.g., worry that if TV-body orbits correspond- 
ing to integrable characteristics have large positive Lya- 
punov exponents, 'regular' flows comprised of such orbits 
would exhibit behaviour resembling chaotic phase mixing 
in smooth potentials! Understanding the bulk properties 
of phase mixing in TV-body systems is important, e.g., in 
light of recent grid code simulations Q indicating that 
chaotic phase mixing may be responsible for 'anomalous 
relaxation' observed in chargcd-particle beams, includin 
the University of Maryland 'five beamlet' experiment 

A complete resolution of these issues will require an 
analysis of 'honest' {i.e., direct summation numerical 
integrations for systems comprised of very large particle 
number TV, which is not yet practical. However, consid- 
erable insight may be derived by considering the simpler 
case of orbits and orbit ensembles evolved in 'frozen- TV' 
systems, i.e., time- independent TV-body systems gener- 
ated by randomly sampling a specified smooth density 
distribution. In particular, by comparing orbits in such 
frozen- TV systems with characteristics with the same ini- 
tial condition evolved in the corresponding smooth po- 
tential, one can quantify the extent to which, as a func- 
tion of TV, TV-body orbits and smooth potential charac- 
teristics actually coincide. Such is the objective of this 
paper. 

Section II focuses on the behaviour of orbits and or- 
bit ensembles in frozen- TV realisations of two simple po- 
tentials, one integrable and the other almost completely 
chaotic. The principal conclusions here are (i) that there 
is a well-defined macroscopic convergence towards the 
continuum limit, both for individual orbits and for or- 
bit ensembles, and (ii) that discreteness effects can be 
extremely well-mimicked by Gaussian white noise in the 
context of a Fokker-Planck or Langevin description. 

Section III considers the more realistic 'thermal equi- 
librium model' [1(1 ]. well known from beam dynamics. 
For appropriate choices of parameter values, this model 
admits yjj large measures of both regular and chaotic 
orbits, so that one encounters a new feature namely 
transitions between regular and chaotic behaviour trig- 
gered by discreteness effects. As for the simpler models 
considered in Section II, one observes clear distinctions 
between regular and chaotic phase mixing, although dis- 
creteness effects, again well modeled by a Fokker-Planck 



description, can be surprisingly important. Even when TV 
is large, individual orbits can exhibit frequent changes in 
behaviour corresponding macroscopically to transitions 
between regularity and chaos; and the scaling implict 
in a Fokker-Planck description suggests that, for chaotic 
orbits, even a total particle number as large as TV ~ 10 9 
may not be large enough to justify a continuum limit. 

Section IV focuses on Lyapunov exponents and the 
meaning of chaos in TV-body systems. The principal con- 
clusion here is that two distinct 'types' of chaos can be 
present in the TV-body problem, characterised by two dif- 
ferent sets of Lyapunov exponents associated with physics 
on different scales. Close encounters between particles 
trigger microchaos, a generic feature of the TV-body prob- 
lem, which leads to large positive Lyapunov exponents 
Xn- If, however, the bulk smooth potential is chaotic, one 
also encounters macrochaos, which is again characterised 
by positive, albeit typically much smaller, Lyapunov ex- 
ponents xs- TV-body realisations of integrable systems 
remain chaotic, even for large TV, in the sense that Xn 
does not decrease towards zero for TV — > oo. Despite this, 
however, microchaos becomes progressively less impor- 
tant macroscopically in that the range of this chaos, i.e., 
the scale on which the microchaos-driven exponential di- 
vergence of nearby orbits terminates, decreases with in- 
creasing TV. 

Section V summarises the principal conclusions and 
speculates on potential implications. 



II. TV-BODY FLOWS AND TV-BODY ORBITS IN 
REGULAR AND CHAOTIC POTENTIALS 

A. Models considered 

The computations described in this Section were per- 
formed for two models which, albeit not representative of 
'real' equilibrium systems, are of significant pedagogical 
value in illustrating the nature of the continuum limit. 
In particular, since one model is integrable and the other 
almost completely chaotic, it is simple to identify sep- 
arately how discreteness effects impact regular versus 
chaotic orbits, an issue that becomes more difficult in 
'realistic' systems which admit a complex coexistence of 
both regular and chaotic orbits. 

Model 1. A constant density triaxial ellipsoid, for which 
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For m < 1, this corresponds to a space-charge potential 
$ sc (r) = $ - 1 (uj 2 a x 2 + u 2 y 2 + uj 2 z 2 ) (3) 
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with $o a constant and frequencies u) a , LUb, and lo c re- 
lated to the axis values a, 6, and c in terms of incomplete 
elliptic integrals. The system was assumed confined by 
an external potential $ ea;t = — 2<£> sc . ft follows that, in 
the continuum limit, each charge evolves in an integrable 
time-independent potential of the form (modulo a con- 
stant) 

W r ) = 2 K* 2 + ^y 2 + u2z2 ) w 

Attention focused primarily on the parameter values a = 
1.95, b — 1.50, and c = 1.05, which, assuming units for 
which Q = 1, implies 13] that w a ss 0.4663, « 0.5508, 
and w c « 0.6753. It follows that the orbital time scale 
t D ~ 2tt/uj ~ 10. 

Because this potential is integrable, one knows that, 
in the continuum limit, only regular phase mixing is pos- 
sible. However, the situation is even more exceptional: 
because of the harmonic character of the potential, i.e., 
the fact that the force is linear, every charge will orbit 
with the same frequencies so that, in the absence of dis- 
creteness effects, there can be no systematic phase mixing 
and no emittance growth in an initially localised clump. 
All emittance growth associated with this potential must 
be attributed to discreteness effects. 

Model 2. Perturbing Model 1 by introducing a spheri- 
cally symmetric, attractive spike of charge near the ori- 
gin, which yields a modified potential 

$c/ la (r) = <S>reg(r) - == (5) 

Vr + i z 

with £ — 10 -3 . Attention here focused on a central 
charge q = 10~ 15 <3 ~ 0.03162, which leads to a potential 
for which, for orbits restricted energetically to m < 1, the 
phase space is almost completely chaotic |l4|. (The bulk 
properties of the potential are insensitive to the precise 
value of I for I < 10 -2 or so; but most of the chaos is 
lost for much larger values of £.) 

Frozen- TV charge density distributions of the form 

1 N 

PN = j^^2 s o{r - ri) (6) 

i=l 

were generated by randomly sampling the smooth density 
distributions p. These correspond to TV-body potentials 

1 N 1 

iV l~i V(r -r,) 2 +e 2 

which incorporate a tiny softening parameter e |l5| . Un- 
less otherwise stated, all Figures in this paper were gen- 
erated from integrations with e = 10~ 5 . 

Orbits were integrated in frozen-iV realisations with 
TV < 10 6 using a variable time step scheme that conserved 
the energy of each charge to at least one part in 10 5 . Esti- 
mates of the largest (finite time) Lyapunov exponent |l6| 



were obtained in the usual way by simultaneously track- 
ing the evolution of a small initial perturbation, periodi- 
cally renormalised at fixed time intervals At. 

The efficacy of phase mixing was tested by generat- 
ing localised clumps of 1600 initial conditions sampling a 
phase space region of size | Ar| ~ | Av| ~ 10 ~ 3 and evolv- 
ing these into the future in both the smooth potential and 
the corresponding TV-body potential (7). The resulting 
orbital data were analysed to compute emittances 



( a = y {rl){vl) - (r a v a ) 2 (a = x,y,z), (8) 
as well as the total 

e=(e x e y e z ) 1 /\ (9) 

The degree to which individual TV-body orbits did, or 
did not, 'look highly irregular' was quantified by com- 
puting the orbital complexity 01 of their Fourier spec- 
tra. This entailed determining for each orbit the quan- 
tities n x , n y , and n z , defined as the minimum number 
of frequencies required to capture a fixed fraction k of 
the power in each direction, and then assigning a total 
complexity 

n(k) — n x + n y + n z . (10) 

In order to obtain a reasonably sharp Fourier spec- 
trum, orbital data were typically recorded at intervals 
St = 0.05, a time corresponding to less than 1% of a 
typical orbital period. 

B. Regular and chaotic phase mixing 

In the continuum limit, initially localised clumps char- 
acterised by the integrable potential (4) exhibit no sys- 
tematic tendency to disperse. Because each orbit ex- 
ecutes harmonic motions with the same three frequen- 
cies, the charges remain close together, returning to their 
original xq, y$, and z$, after periods t x , t v , an t z . Dis- 
creteness effects break this exact periodicity and trigger a 
systematic spread. This is illustrated in Figure 1, which 
exhibits the x and y coordinates of the same 1600 orbit 
ensemble with E — 1.0 at five different times, allowing 
for frozen- TV backgrounds with N = 10 3 5 and N = 10 5 . 
For N — 10 3 5 , this dispersal is comparatively rapid, the 
charges having spread to sample the entire allowable con- 
figuration space within t = 128, a time corresponding to 
only 10 orbital periods or so. For the larger system with 
N = 10 5 , the dispersal is considerably slower, requiring 
a time t ~ 512, roughly four times larger, to achieve a 
comparable spread. 

The situation is very different for the potential (5), for 
which, even in the continuum limit, the particle phase 
space is almost completely chaotic. In this case, one 
observes exponentially fast chaotic phase mixing in the 
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smooth potential, and allowing for discreteness effects 
only accelerates the process. This is evident from Figure 
2, which exhibits the x and y coordinates for a 1600 orbit 
clump, again with E = 1.0, evolved in frozen- N realisa- 
tions with N = 10 3 ' 5 , N = 10 4 - 5 , and N = 10 5 5 , as well 
as (in the right hand column) the smooth potential. Even 
for the smooth potential, a time t ~ 128 is sufficient for 
particles to sample most of the energetically accessible 
phase space. 

The visual impression that the chaotic clump disperses 
far more rapidly can be quantified by computing the 
omittance e as a function of time. The left hand column 
of Figure 3 exhibits e(t) for the same ensemble of initial 
conditions used to generate Figure 1, now allowing for 
frozen-TV backgrounds with N = 10 3 , 10 3 ' 5 , 10 4 , 10 4 5 
and and 10 5 . For the smallest value of N it is not com- 
pletely clear whether e grows exponentially or as a power 
law in time. However, for N > 10 3,5 , the growth is dis- 
tinctly subcxponential. Indeed, the data for N > 10 3 5 
are well fit by an emittance growth law 

e^it/to) 1 ' 2 (11) 

where 

t G oc Nt D . (12) 

Extrapolating to the limit N — > oo yields the expected 
result that there can be no systematic emittance growth. 

The left hand column of Figure 4 exhibits analogous 
results for the initial conditions used to generate Fig- 
ure 2, now plotted on a logarithmic scale, allowing for 
N = 10 2 ' 5 , 10 3 - 5 , 10 4 - 5 , 10 5 - 5 , and, in the bottom panel, 
the smooth potential. It is evident that, for the smooth 
potential, In e exhibits a roughly linear growth during the 
interval (say) 10 < t < 100, corresponding to an expo- 
nential growth in emittance. This is hardly surprising. 
The fact that individual orbits in the clump are chaotic 
implies that they should diverge exponentially so that, at 
least for small e(0), one would expect e to grow exponen- 
tially at a rate comparable to the value of a typical (finite 
time) Lyapunov exponent xs for the smooth potential. 
For this ensemble, the mean exponent for the interval 
< t < 256 assumed the value (xs) = 0.056, which 
corresponds to the slope of the dashed line in panel (i). 

Allowing for discreteness effects clearly accelerates the 
rate of chaotic phase mixing. For the two smaller values 
of N, 10 2 - 5 and 10 3 5 , the growth again appears expo- 
nential, albeit at a larger rate; but for the systems with 
N = 10 4 ' 5 and 10 5 5 the evolution is clearly more com- 
plex. Indeed, a careful examination of the data for these 
two cases suggests strongly that the evolution can be de- 
composed into two largely distinct exponential phases, 
the former characterised by a growth rate 3> Xs ar >d the 
latter by a slower rate ~ xs- This is consistent with the 
analysis to be presented later in Section IV, which indi- 
cates that two sorts of chaos can act in A-body systems, 



large scale macrochaos characterised by a Lyapunov ex- 
ponent xs and shorter range microchaos characterised by 
a Lyapunov exponent xn S> Xs- (For very small A, the 
microchaos also acts on macroscopic scales, thus over- 
whelming any observational effects associated with the 
much weaker macrochaos: hence the (near-) absence of 
the second exponential phase in panels (a) and (c)!) The 
data summarised in Figure 4 are consistent with a second 
exponential phase with 

ecx AT-^exptxs*), (13) 

the form of which will be motivated in Section lid. 

The results derived here for Model 2 which, in the 
continuum limit, corresponds to a chaotic potential, are 
likely generic for bulk density distributions correspond- 
ing to a chaotic potential. However, the results for Model 
1 are special in that there is no systematic emittance 
growth in the continuum limit. If, as one would expect 
in a 'real' system, the bulk potential exhibits at least 
some anharmonicities, regular phase mixing will trigger 
linear emittance growth even in the continuum limit. In 
this case, allowing for discreteness effects will again ac- 
celerate the growth of emittance, but the exact form of 
this enhanced growth can be more complex, even though 
it will again be subexponential. Section III will exhibit 
additional examples of how discreteness can accelerate 
emittance growth for both regular and chaotic clumps. 

C. Individual Orbits and the Continuum Limit 

The preceding indicates that, as A increases, orbit en- 
sembles evolved in frozen-TV backgrounds more closely 
resemble orbit ensembles with the same set of initial 
conditions evolved in the corresponding smooth poten- 
tial. This does not, however, necessarily imply that indi- 
vidual orbits also converge towards characteristics in the 
smooth potential. To what extent, then, is it true that, as 
A increases, individual trajectories come to more closely 
resemble smooth potential characteristics? 

The most obvious - and compelling - check is visual: 
do frozen-A orbits 'look like' smooth potential charac- 
teristics when A becomes sufficiently large? As a sim- 
ple, and extreme, example, consider a constant density 
spherical system without a central spike where, in the 
continuum limit, $ reg reduces to an isotropic harmonic 
oscillator potential; and select an initial condition which, 
in the smooth potential, corresponds to a circular orbit. 
Results derived from integrations of such an initial con- 
dition in different frozen- A systems are exhibited in Fig- 
ure 5, which shows representative frozen- A'' orbits gener- 
ated for particle number varying between A = 10 2 5 and 
N = 10 5 ' 5 , along with the smooth potential orbit. For 
the four smallest values of A, there is no obvious hint 
that the orbit 'should' be circular or that there 'should' 
be a net sense of circulation, although there is a crude 
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visual sense that, as A increases, the orbit becomes 'less 
tangled.' However, for TV = 10 4 ' 5 one starts to discern a 
clear sense of net circulation, for A = 10 5 the orbit has 
clearly become centrophobic (thus suggesting that angu- 
lar momentum is at least approximately conserved), and 
for A = 10 5 ' 5 the orbit arguably resembles a 'distorted' 
or precessing circular orbit. 

As is illustrated in Figure 6, the visual impression that 
the orbit is becoming more nearly circular can be cor- 
roborated by constructing the Fourier spectra of the or- 
bital data. For the three smallest values of A, the quan- 
tity \x(uj)\ is obviously broad band, although there is a 
peak at or near the circular frequency associated with 
the smooth potential orbit. For A = 10 40 and 10 4 5 the 
peak becomes appreciably sharper, and for A = 10 5 
and 10 5 5 one sees only slight irregularities in the spectra 
which translate into the visual appearance of precession. 

The conclusion is obvious: as A increases, frozen-TV 
orbits come to more closely approximate the smooth or- 
bit, both visually and in terms of their power spectrum. 
Analogous results obtain for more generic initial condi- 
tions evolved in this and other integrable potentials. 

Convergence of orbits in terms of their Fourier spec- 
tra is important in justifying straightforward applications 
of nonlinear dynamics to many-body systems interact- 
ing via long range forces. Many physical phenomena 
in many-body systems, including accelerator modes 181, 
modulational diffusion |l9| , and resonant relaxation 2^ , 
are attributed to resonant couplings between, e.g., the 
natural frequencies of individual regular orbits and the 
frequency or frequencies associated with some perturba- 
tion. However, such applications can only be justified if 
the 'real' iV-body orbits have frequencies that adequately 
approximate the frequencies associated with characteris- 
tics in the smooth potential. 

The degree of irregularity exhibited by individual or- 
bits can be quantified by computing Fourier complexity, 
as defined in Eq. (10). The results of one such investi- 
gation are summarised by the curves with diamonds in 
Figure 7, which exhibit n(0.95), the mean number of fre- 
quencies required to capture 95% of the total power, for 
collections of 100 initial conditions evolved in frozen- A'' 
backgrounds with different A. (The triangles will be dis- 
cussed in the following subsection.) In each case, the ini- 
tial conditions were integrated for a time T = 128, with 
orbital data recorded at fixed intervals St — 0.05. The 
data were then Fourier analysed using an FFT solver to 
translate a set of 2j points into a set of j Fourier am- 
plitudes. The upper panel was generated for orbits in 
the integrable Model 1, the lower panel for the strongly 
chaotic Model 2. In each panel, the solid curve repre- 
sents the mean complexity computed in the unperturbed 
smooth potential. Error bars were computed by dividing 
the 100 initial conditions in half and analysing each half 
separately. 



ing function of A which converges towards the continuum 
value for A — > oo. For smaller values of A, the regular 
and chaotic ensembles have comparable complexities, al- 
though the regular ensemble is slightly less (~ 20%) com- 
plex. However, for larger values of A there are clear dis- 
tinctions between the regular and chaotic ensembles, the 
chaotic ensembles for A > 10 5 being nearly twice as com- 
plex as the regular ensembles. Indeed, for other choice 
of regular models the complexity can be even lower: The 
fact that the smooth potential complexity in panel (a) is 
as large as it is reflects the fact that the regular orbits 
used to generate this Figure were relatively complex 'box' 
orbits, with the topology of Lissajous figures, which, even 
in the continuum limit, require two or three frequencies in 
each direction to capture as much as 95% of the power. If 
instead n(0.95) is computed for an ensemble of initial con- 
ditions corresponding to circular orbits, the complexity 
converges towards a continuum limit with ra(0.95) = 3. 



D. Modeling iV-body orbits and flows by Gaussian 
white noise 



Conventional wisdom holds that discreteness effects 
can be idealised as friction and Gaussian (nearly) white 
noise in the context of a Fokker-Planck description fj. 
Taken literally, this suggests that individual A-body or- 
bits can be well-mimicked by Langevin simulations. How- 
ever, it is not completely clear to what extent this is really 
true. The original derivation of the Fokker-Planck equa- 
tion (and most if not all of its tests) restricted attention 
to the statistical properties of orbit ensembles or distri- 
bution functions over comparatively long time scales, as- 
suming implicitly that the bulk potential in which the 
particle evolves is nonchaotic. Can the friction/noise 
paradigm describe correctly short time behaviour and/or 
the behaviour of individual orbits, especially in a chaotic 
potential? 

Granted the validity of a Fokker-Planck description, 
a simple rule connects A to the strength of the friction 
and noise. Assuming that the noise is characterised by a 
temperature per unit mass comparable to the magni- 
tude of the particle energy, the coefficient of dynamical 
friction r\ defines an energy relaxation time tji = r/^ 1 . 
However, an evaluation of the Fokker-Planck coefficients 
in a binary encounter approximation leads to the predic- 
tion that tR oc (A/lnA)<£>, with In A the Coulomb log- 
arithm. Given the assumption of a nonneutral plasma, 
the treatment of A must necessarily be somewhat heuris- 
tic |2lj. However, there is a general agreement that A 
should scale as some power of A, so that tn, and hence 
the diffusion constant D, should satisfy 



t R oc D oc 77 oc 



In A 
~A~' 



(14) 



It is evident that, in both cases, n(0.95) is a decreas- The obvious question, then, is whether frozen-A Simula- 
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tions with specified N can be well-mimicked by Langevin 
simulations with r\ cx (In N/N). 

Two practical issues arise in testing this prediction. 
The first is that, because of the limited range of N that 
can be explored, it is impractical to test the subdominant 
In AT dependence: for iV < 10 3 or so, the very notion of 
a bulk potential fails; for N > 10 6 computations become 
prohibitively expensive. One must instead restrict at- 
tention to testing the simpler scaling relation 77 oc N , 
i.e., 

In r) = p — In N, (15) 

for some constant p. 

The second point is more serious. The usual Langevin 
equation reads [4| 

ff^ v civ 

^ = _V *-T^+J' o , (a = x,y,z), (16) 

where rjdr a / dt represents a dynamical friction. F a repre- 
sents Gaussian white noise, which is characterised com- 
pletely by its first two moments: 

(F a (t))=0, {a,b = x,y,z) 

and 

(F a (h)F b (t 2 )) = 2r 1 G6 ab 5 D (t 1 - t 2 ), (17) 

with D = 2r]Q the diffusion constant entering into a 
Fokker-Planck description. By choosing 9 to equal the 
initial energy one can ensure that the average energy of 
the orbits remains unchanged. 

Such an equation is clearly unsatisfactory here. En- 
ergy is conserved absolutely for frozen- iV orbits, so that 
one must also impose energy conservation on any scheme 
which aims to mimic its effects. (For very small 77, energy 
remains almost conserved for very long times. However, 
comparatively small ./V should correspond to relatively 
large 77, which implies large changes in energy and, as 
such, significant changes in the phase space regions ac- 
cessible to the noisy orbit.) For this reason, the noisy 
integrations described here were performed using a mod- 
ified energy-conserving noise [z| • 

This entailed (1) eliminating the dynamical friction al- 
together, (2) again imparting random kicks as in Eq. (17), 
but (3) renormalising the modified velocity at each time 
step by an overall factor, i.e., v(t + St) — > av(t + St), 
with a so chosen that E(t + St) = E(t). Modulo this 
complication, the noise was integrated using a standard 
algorithm |23j based on a fourth order Runge-Kutta in- 
tegration scheme with fixed time step St. The integra- 
tions were performed for St = 2 x 10~ 4 , it having been 
confirmed that the statistical effects of the noise were 
insensitive to the precise value of St for St < 10 -3 . 

At the level of orbit ensembles, as probed by the emit- 
tance and other bulk moments, the results of frozen- 
N simulations are in fact extremely well-mimicked by 



Langevin simulations, at least for comparatively large N. 
The degree to which this is true can be inferred by con- 
trasting the right and left hand columns of Figures 3 and 
4. As discussed already, the left hand columns of Figures 
3 and 4 exhibit, respectively, time-dependent emittances 
for Models 1 and 2, allowing for several different values 
of N. The right hand panels exhibit data generated from 
Langevin integrations of the same initial conditions, al- 
lowing for amplitudes 7/ satisfying Eq. (15) with p = 0.5, 
so that, e.g., N = 10 5 5 corresponds to 7/ = 10~ 50 . For 
the smallest values of N (and hence the largest values of 
77) - corresponding to panels (a) and (b) in Figure 3 and 
panels (a) - (d) in Figure 4, the agreement is not all that 
good. However, for larger particle number - N > 10 3 5 
for the regular system and N > 10 4 5 for the chaotic sys- 
tem -, the agreement is obviously quite good. 

The bottom right hand panel in Fig. 4 was generated 
for orbits evolved with a considerably smaller value of 
77, namely rj = 10~ 7 ' 5 , this corresponding to the largest 
noise amplitude that does not alter appreciably the emit- 
tance growth observed in the smooth potential. To the 
extent that the scaling of Eq. (15) is in fact correct for 
p « 0.5, the fact that larger values of 77 have an apprecia- 
ble effect on emittance growth implies that, even over an 
interval as short as t = 128, corresponding to ~ 10 orbital 
time scales tjj, one requires A~ > 10 7 to justify a contin- 
uum limit! Even though the collisional relaxation time 
scale tn tx (N/\nN)tjj 3> tp, discreteness effects can be 
important in a system with N ~ 10 6 ' 5 ott, a time scale as 
short as ~ lOtp. 

As is evident from Figure 7, this agreement also ex- 
tends to the level of individual orbits. As described 
already, the diamond curves in this Figure were de- 
rived from frozen-Af integrations. The other curves, 
constituted of triangles, were generated from exactly 
the same initial conditions, now integrated, however, in 
the smooth potential while being subjected to energy- 
conserving Gaussian noise with 8 = E and variable rj. To 
the extent that conventional Fokker-Planck theory is cor- 
rect, one would anticipate a correspondence between N 
and 77 of the form given by Eq. (15). The noisy points in 
Figure 7 were in fact identified with the frozen- A" points 
assuming the validity of the scaling (15) with p = 0.6. 
The obvious fact, then, is that, given this identifica- 
tion, the curves n(N) and 77,(77) rather nearly coincide. 
Even at the level of the complexity of individual orbits, 
frozen-N orbits can be well-mimicked by noisy orbits with 
In 77 + In A^ = const. 

Granted that discreteness effects can be mimicked by 
Gaussian white noise, the scaling relations (11) and (13) 
are easily understood. At least for a harmonic potential, 
it is straightforward to derive analytic solutions to the 
Langevin equation (16) for moments like (x 2 ) or (xv x ) 0. 
Alternatively, it is easily seen that the Fokker-Planck 
equation associated with Eq. (16) implies that the clump 
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emittance satisfies 

^ = 2 V + 2 V ((x 2 )(v 2 )-(xv) 2 ). (18) 

Assuming, however, that the initial emittance is ex- 
tremely small, at early times one can approximate that 
(x 2 ) (v 2 ) w (xv) 2 ; and, to the extent that the growth time 
is long compared with the characteristic crossing time, 
one can average over oscillations to set (x 2 ) — E/w 2 , 
with E the initial energy. It then follows that, for early 
times, 



/2J5B/// ^ 1 2 



(19) 



Combining Eq. (19) and the analogous formulae for e y 
and e z with the scaling relation r\ oc 1/N leads immedi- 
ately to Eq. (11). The same diffusive t 1 / 2 behaviour also 
arises for Sr rms and Sv rms - 

A somewhat more heuristic argument can account for 
the scaling (13) associated with a chaotic potential. If 
the initial emittance e(0) — 0, it is clear that, in the ab- 
sence of discreteness effects, e(t) would continue to vanish 
identically: two smooth integrations of the same initial 
condition will not yield divergent orbits, even if the orbits 
are chaotic. However, discreteness effects act to 'kick' two 
nearly coincident orbits apart, at which point they will 
tend to diverge at a rate set by the Lyapunov exponent 
Xs associated with the bulk potential. Assuming, how- 
ever, that the 'kicks' are random, their effects will scale 
as 77 1 / 2 rather than 77; but combining this with Eq. (15) 
implies that 



$v rms oc Sv rrns oc N-^ 2 eMxst). 
Eq. (13) follows since e scales as Sr rms and 5v r 



(20) 



III. THE THERMAL EQUILIBRIUM MODEL 
A. Defining the model 

Consider now a more realistic example, the thermal 
equilibrium model jlOj . which, in the continuum limit, 
admits large measures of both regular and chaotic or- 
bits 0|- This model allows for a collection of N identi- 
cal charged particles, interacting electrostatically, that is 
constrained by linear restoring forces to manifest triax- 
ial symmetry, the focusing forces in different orthogonal 
directions being characterised in general by unequal fre- 
quencies. Individual particles thus have energy 



E = -m» 2 + — to(oj-x) 2 + 1 



(21) 



where x and v denote particle position and velocity, m 
and q denote the mass and charge, u) — (u> x ,uj y ,uj z ) rep- 
resents the three frequencies associated with the focusing 
force, and </>(x) is the collective space-charge potential. 



The additional assumption is that the particles can be 
characterised by a one-particle distribution function of 
the Maxwell-Boltzmann form, / oc exp(— E/ksT), with 
fcflT the temperature. This implies a bulk number den- 
sity satisfying 



n(x) = n(0) exp 



-4)m<(ui-x.) 2 — q4>(x) 



(22) 



where 0(x) is defined implicitly as a function of n via the 
relations (in mks units) 

V 2 0(x) = -— n(x), (f)(0) = V <f>(0) = 0. (23) 

Following, e.g., ^lj> the problem can be cast into di- 
mensionless form by expressing length and frequency in 
units of the Debye length and plasma frequency, i.e., 



A i 



e k B T 



n(0)q 2 



n(0)q 2 ' *" £o m 

and by introducing a dimensionless potential 



$(x) = 



k B T 



(24) 



(25) 



With appropriate rescaling, the net result is a density 
distribution of the form 



ra(x) = exp 



jft 2 i? 2 (x)-$(x) 



(26) 



where 



V 2 $(x) = -n(x), $(0) = V$(0) = 0. (27) 

Here il 2 = (uo y /Ld p ) 2 and R 2 = (x/a) 2 + y 2 + (z/c) 2 , in 
terms of scale lengths a and c satisfying a — u) y /u> x and 
c = u)y/u) z . The minimum permissible focusing strength 
corresponds to 



n = n u = i/V(i/a 2 ) + i + (i/c 2 ). 



(28) 



The experiments described here were performed assum- 
ing parameter values a 2 = 0.5, c 2 = 1.5, and Q = 
1.0001/\/3, for which a typical orbital time scale to ~ 20. 

These parameters represent a beam that is moder- 
ately, but not strongly, dependent on space charge. Con- 
sider, for example, a proton bunch with 1 /Ltm root mean 
squared normalised emittance spanning 3 cm full 'ra- 
dius'. If the bunch is described by the thermal equilib- 
rium model, the Debye length is ~ 2 mm and the bunch 
population is the ~ 3 x 10 9 protons, this corresponding 
to a bunch charge ~ 0.5 nC. 

In general it does not appear possible to solve eqs. (26) 
and (27) analytically. This makes both the generation of 
iV-body realisations of the density and the computation 
of orbits in the smooth potential much more difficult. 
However, these difficulties can be, and were, resolved 
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using numerical techniques described in In princi- 

ple, the accelerations for the N-body thermal equilibrium 
model should be given by 



-n 2 R 2 



1 TV 

4^TVl 



jv 

E 



(29) 



with TV the number of frozen particles and TV m satisfying 



N m = 



d 3 rexp 



1 



-n 2 R 2 - $(r) 



(30) 



In practice, however, one cannot perform this integral, 
even numerically, since $ was only evaluated on a fi- 
nite grid. For this reason, the integral was first solved 
with limits coinciding with the grid boundaries, but then 
renormalised by a small constant factor so that plots of 
the potential and density in the smooth and frozen-TV 
configurations overlapped perfectly. 



B. Regular and chaotic phase mixing 

The top four left hand panels of Fig. (8) exhibit emit- 
tance growth for an initially localised orbit ensemble 
evolved in frozen-TV realisations of the thermal equilib- 
rium model, selected with energy sufficiently small that 
the constant energy hypersurface in the smooth poten- 
tial is completely regular. Since the size of the accessible 
phase space is roughly ten times larger than was the case 
for the oscillator models, the initial conditions sampled a 
region ten times larger, |Ar| ~ |Au| ~ 10~ 2 . It is evident 
that, as for the integrable oscillator model considered in 
Section II, the emittance growth is power law rather than 
exponential; and, at least for TV — 10 5 5 and TV = 10 6 , 
it is well fitted by a growth law e cx t 1 ^ 2 . 

The bottom left panel exhibits emittance growth for 
the same orbit ensemble evolved in the smooth poten- 
tial. Here the evolution is clearly linear, rather than 
square root, the expected behaviour for smooth orbits 
in generic integrable potentials where nearby initial con- 
ditions correspond to slightly different orbital frequen- 
cies. (Recall that, for the oscillator Model 1, there is 
zero emittance growth in the continuum limit.) That e 
grows much faster for the frozen- TV model with TV = 10 60 
than for orbits in the smooth potential indicates clearly 
that, for the thermal equilibrium model, TV = 10 6 is not 
a sufficiently large particle number to justify a continuum 
approximation, even over an interval as short as t = 512. 

The top four right hand panels demonstrate that dis- 
creteness effects can again be well-mimicked by energy- 
conserving Gaussian white noise with 6 = E and TV 
and r\ related as in Eq. (15) although, in this case, the 
best fit value p w 1.5, rather than p» 0.5. The bottom 
right panel exhibits the emittance growth for an ensem- 
ble evolved with r\ = 10~ 65 , the largest noise amplitude 
that does not significantly alter emittance growth in the 



smooth potential. Presuming that the scaling (15) can 
again be extended to larger TV and smaller rj, one infers 
that, for the purpose of predicting emittance growth in 
this regular ensemble, the smallest value of TV for which 
the continuum limit can be justified is TV ~ 10 8 . 

Figure 9 exhibits analogous data for a higher energy en- 
semble which, in the continuum limit, corresponds com- 
pletely to chaotic orbits. It is evident that, as for the 
chaotic model in Section II, the evolution is exponential 
overall, rather than power law; and that discreteness ef- 
fects again have an important effect. Also evident from 
a comparison of left and right hand panels is that, as for 
regular orbits, discreteness effects can be well-mimicked 
by noise with TV oc 1/rj and p sa 1.5. Most striking, how- 
ever, is the fact that, in this case, even much weaker noise 
can accelerate emittance growth appreciably. For this 
chaotic ensemble, discreteness effects must correspond to 
a noise amplitude satisfying rj < 10~ 80 or so before a con- 
tinuum limit can be justified. Chaotic orbits are far more 
susceptible to low amplitude noise than are regular orbits. 
Presuming again that the scaling relation (15) holds, this 
implies that for the case of chaotic orbit ensembles, the 
continuum limit cannot be justified for TV < 10 9 ' 5 . 

This, coincidentally, is roughly the number of parti- 
cles in the equilibrium proton bunch described in Sec. 
Ill A. Accordingly, in studying the dynamics of beams 
with moderate space charge, one may not be able to as- 
sume the validity of the continuum limit with complete 
confidence, even for a system in equilibrium. The situ- 
ation may be even more problematic for a beam that is 
significantly out of equilibrium, since the resulting time- 
dependent potential would be expected to generate a 
larger population of chaotic orbits . 



Transitions between regular and chaotic 
behaviour 



At very low energies, where the total potential is nearly 
harmonic, all smooth potential orbits are regular, so that 
discreteness effects can only act to deflect frozen- TV or- 
bits from one regular trajectory to another. However, at 
higher energies the smooth potential admits a complex 
coexistence of regular and chaotic orbits. This implies 
the possibility that discreteness effects can deflect frozen- 
TV orbits from regular to chaotic characteristics and vice 
versa _12|. Of obvious interest then is how fast, as a 
function of TV, such transitions occur. 

For example, an accelerator designer relying on the 
Vlasov equation and an analysis based on a smooth, 
macroscopic potential would neglect these microscopic 
transitions and their corresponding impact on chaotic 
mixing. Thus the physics of collective relaxation and 
global emittance growth would be improperly modeled 
and the results, at least in principle, would be suspect. 
The consequences of such an omission depend on the 
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problem at hand, but one might expect them to be es- 
pecially severe for nonequilibrium beams where chaotic 
dynamics is likely to be more prevalent [25|. 

If one selects a localised ensemble of initial conditions 
corresponding to regular orbits in the smooth potential 
and integrates these initial conditions into the future, 
discreteness effects will, if sufficiently strong, eventually 
trigger transitions from regularity to chaos. That such 
transitions actually occur can be determined by a visual 
inspection of individual frozen- ./V orbits which can be ob- 
served to become abruptly 'more irregular' in appear- 
ance. If, moreover, large numbers of transitions occur 
over very short times, this can make the emittance asso- 
ciated with an initially localised ensemble, which ought 
to grow as a power law in time, exhibit instead a more 
rapid, roughly exponential, increase. 

However, an accurate determination of the relative 
fraction of the orbits which are still regular requires a 
more careful analysis. This was done by recording the 
phase space coordinates of individual frozen- TV orbits at 
various times t > 0, and evolving these into the future 
in the smooth potential to determine whether the result- 
ing smooth potential characteristics were still regular or 
whether instead they had become chaotic. The most 
straightforward fashion in which to determine whether 
the orbits are chaotic would have been to compute an 
estimate of the largest (finite time) Lyapunov exponent. 
Given, however, that the potential cannot be expressed 
in terms of elementary functions, this would have proven 
extremely expensive computationally. For that reason, 
distinctions between regularity and chaos were based in- 
stead on the computed complexities of the characteris- 
tics. As discussed elsewhere (e.g., 01 i E3> an d refer- 
ences cited therein), such a criterion typically coincides 
almost exactly with more conventional criteria based on 
Lyapunov exponents. 

Presuming that the system is ergodic and that dis- 
creteness effects are sufficiently strong that they can in 
principle convert any orbit from regular to chaotic, and 
vice versa, it would seem clear what such an analysis 
ought to reveal. (1) At sufficiently late times, indepen- 
dent of TV the relative fraction of chaotic orbits generated 
from any initial ensemble should (to within statistical un- 
certainties) coincide with the relative measure of chaotic 
orbits on the constant energy hypersurface, i.e., to the 
relative volume of the chaotic portions of the constant 
energy hypersurface. (2) Assuming, however, that dis- 
creteness effects are more important for smaller TV, the 
time required to converge towards this asymptotic value 
should be an increasing function of TV. As TV increases, 
transitions should become more rare. 

As illustrated in Fig. 10, which summarises computa- 
tions with particle number between 10 4 5 and 10 6 , this 
expectation was in fact confirmed. For frozen-TV sys- 
tems with number as small as TV = 10 4 ' 5 , nearly 40% of 
the orbits in an initially regular ensemble had become 



chaotic within a time t = 64, a time corresponding to 
only ~ 3<d, and the relative measure / of chaotic orbits 
appears to have asymptoted towards a time-independent 
value by t = 128. The relative measure of chaotic orbits 
for computations with TV = 10 5 ' grew more slowly in 
time; but, by t — 512, the relative measure had again 
approached a value / ~ 0.4. For larger values of TV, f re- 
mains a monotonically increasing function of time, but 
transitions are sufficiently rare that one does not ap- 
proach an equilibrium population within a time as short 
as t = 512 ~ 25t D . 

IV. LYAPUNOV EXPONENTS FOR 
MICROCHAOS AND MACROCHAOS 

A. Ordinary Lyapunov exponents 

As TV increases, frozen- TV orbits come to more closely 
resemble smooth potential characteristics generated from 
the same initial condition, both visually and in terms of 
their Fourier spectra. One might therefore expect that, 
at least for a regular, i.e., nonchaotic, smooth potential, 
the value of the largest Lyapunov exponent xn should 
decrease with increasing TV and converge towards zero 
for TV — + oo. Such, however, is not the case. Rather, as 
is also true for gravitationally interacting systems of par- 
ticles 22], for both regular and chaotic orbits the value 
of xn is comparatively insensitive to TV; and there are 
even indications that xn might increase with increasing 
TV. 

Figure 11 exhibits the value of the largest Lyapunov ex- 
ponent Xjv as a function of TV for a single initial condition 
evolved in frozen- TV integrations with softening parame- 
ters varying between e = 10 -5 and e = 10 _1 . Figure 12 
exhibits the same data, now plotting xn as a function of 
e for different values of TV. In both Figures, each point 
was generated by integrating the same initial condition 
used to generate Figures 5 and 6 for a total time t = 256 
in 10 different frozen- TV realisations of Model 1. In the 
continuum limit this initial condition corresponds to an 
integrable circular orbit with vanishing Lyapunov expo- 
nents; and, as was evident from Figure 5, the frozen-TV 
orbits for larger TV look much more regular in appear- 
ance than do the orbits with smaller TV. Despite this, 
however, at least for the smallest values of e, xn does 
not decrease with increasing TV. As probed by xn, orbits 
with TV = 10 2 ' 5 and 10 5 are comparably chaotic! 

However, for larger values of e, xn does decrease with 
increasing TV. That this should be the case is easily un- 
derstood. Because the bulk potential is integrable, the 
chaos must at some level be associated with close en- 
counters between nearby charges; but introducing a soft- 
ening parameter de facto 'turns off' encounters on scales 
shorter than e. If the charge density is sufficiently large 
that encounters with separation < e become common, 
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the simulation will have artificially, and erroneously, re- 
duced this source of chaos, yielding a smaller xn- 

That the value of xn f° r nearly unsoftened frozen- 
N integrations is insensitive to whether the smooth po- 
tential is regular or chaotic is illustrated in Figure 13, 
which exhibits xjv as a function of N for integrations 
with e = 10~ 5 for the same initial condition integrated 
in both Models 1 and 2. That Figure also emphasises 
another important point, namely that xn is typically 
larger than any Lyapunov exponent xs associated with 
motion in the smooth potential by an order of magnitude 
or more. For this particular initial condition, xn ~ 0.82 
whereas xs ~ 0.056. 

That xn should be insensitive to the choice of N, at 
least for unsoftened simulations, might seem rather cu- 
rious. However, it is not hard to understand why this 
might be so for a 1/r 2 force. Given that the microchaos 
disappears completely in the continuum limit, it would 
seem clear that it must be associated with a sequence of 
'random' interactions between a 'test' charge and a col- 
lection of 'field' charges. However, this suggests that the 
Lyapunov time = x^ 1 associated with the growth of 
a small initial perturbation can be estimated by consid- 
ering tidal effects associated with a pair of charges sepa- 
rated by a distance r comparable to (some fixed fraction 
of) the typical interparticle separation r sep . This tidal 
acceleration scales as 



d 2 5r _ N q 5r 

— = (<5r-V)aoc — 5r = ^, 



dt 2 



sep 



(31) 



with q the magnitude of an individual charge. Given, 



however, that r sep ~ n 1 ^ 



N l ^R sys , with R sys the 



size of the system and n a characteristic number density, 
and that q — Q/N, with Q the total charge, it follows 
that the time scale £*, and hence xn, should be indepen- 
dent of N. As N increases, the sizes of the individual 
charges and the cube of the typical interparticle separa- 
tion both decrease as TV -1 , so that the ratio is indepen- 
dent of particle number. 

A more careful argument |2(| actually allows one to 
prove analytically that, in the absence of softening, xn 
cannot decrease towards zero with increasing N. For 
simple geometries, an analytic expression for the average 
value of the stability matrix entering into the definition 
of xn can be formulated in terms of a 3iV-dimensional 
integral. This integral cannot be evaluated analytically, 
but one can derive rigorous bounds which ensure that 
the largest eigenvalue remains positive even for N — ► oo; 
and, even more strikingly, Monte Carlo evaluations of the 
integrals suggests that xn should be a slowly increasing 
function of N. In other words, viewed in terms of xn, 
orbits become more chaotic as N increases, even though, 
for the case of a regular potential, they become more 
regular in appearance! [13] 

The obvious inference is that N-body Lyapunov expo- 
nents xn do not provide a useful characterisation of the 



degree of chaos associated with an orbit, at least when 
that orbit is viewed macro scopically. 



B. Microchaos and macrochaos in the iV-body 
problem 

That frozen- N orbits have large positive Lyapunov ex- 
ponents xn, even for the case of an integrable potential, 
but that distinctions between regular and chaotic poten- 
tials are clearly manifested in the phase mixing of ini- 
tially localised clumps would suggest that there are two 
different, and comparatively distinct, potential sources of 
chaos in the iV-body problem. 

On the one hand, one would expect microchaos, trig- 
gered by close encounters between individual charges, 
which is manifested only on very short scales, comparable 
to, or perhaps somewhat larger than, the typical inter- 
particle spacing. This source of chaos should be generic 
to the iV-body problem, irrespective of the form of the 
bulk potential, generating randomness qualitatively sim- 
ilar to what arises in pinball. On the other hand, there 
is also the possibility of larger scale macrochaos, which 
would be expected if and only if the bulk potential admits 
global stochasticity. 

If these expectations are in fact correct, two nearby 
initial conditions evolved in a frozen- N realisation of any 
potential should diverge exponentially at a rate ~ xn un- 
til their separation becomes somewhat larger than a typ- 
ical interparticle spacing, at which point the microchaos 
would 'turn off.' If the bulk potential is regular, no other 
source of chaos could act and the two orbits would con- 
tinue to diverge as a more modest power law. If, however, 
the bulk potential is chaotic, macrochaos would still act, 
resulting in a continued exponential divergence, albeit at 
a rate ~ xs typically much smaller than xn- In this case, 
exponential divergence should be replaced by a power law 
divergence only once the separation has become macro- 
scopic. 

This three-stage evolution for chaotic orbits is clearly 
illustrated in Figure 14, which exhibits data for frozen- iV 
simulations with N = 10 5 , 10 5 5 , and 10 6 . The three 
curves in that figure were each generated by selecting 
50 fiducial initial conditions in a phase space region of 
size 5 x 10~ 4 and 50 perturbed initial conditions that 
were subjected to a displacement Sx = 10~ 5 , evolving 
each initial condition into the future and computing as a 
function of time the phase space separation [29j 



5Z= ^\Sr\ 2 + |<5v| s 



(32) 



That SZ experiences two distinct stages of exponential 
evolution, at very different rates is especially evident in 
the curve with N = 10 6 . The solid lines accompanying 
that curve have slopes 0.82 and 0.056, corresponding, re- 
spectively, to the mean values of the TV-body xn and the 
smooth potential xs for those initial conditions. 
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As TV increases, the initial exponential phase termi- 
nates for smaller values of SZ. The microchaos responsi- 
ble for this first phase will 'turn off' when \6v\ becomes 
large compared with a typical interparticle spacing, but 
that interparticle spacing scales as TV -1 / 3 . 

Figure 15 shows the analogue of Figure 14, now gen- 
erated for regular initial conditions in the thermal equi- 
librium model. Once again there is an initial exponential 
growth at a rate comparable to xn, but in this case there 
is no evidence of a second exponential phase. Rather, the 
initial exponential phase is followed immediately by an 
interval of power law divergence. 

One other point, not obvious from these Figures, is 
that the scaling with TV observed for the final power law 
phase is different for the regular and chaotic systems. 
For regular orbits, the phase space separation SZ ', like 
\Sr\ and |£v|, satisfies a linear growth law 

SZ(t) cx (t/t G ), (33) 

with 

t G oc N- l 'H D . (34) 

For chaotic orbits, one finds that SZ again grows linearly 
in time, but that the growth time 

t G oc (l/\nN)t D . (35) 

C. Alternative interpretations of Lyapunov 
exponents 

The standard definition of Lyapunov exponents implies 
that they probe the average rate of divergence for two 
nearby chaotic orbits in a single system. However, for 
the TV-body problem, Lyapunov exponents also quantify 
two other effects which, as a practical matter, are of equal 
importance: 

1. Lyapunov exponents probe the rate at which orbits 
generated from the same initial condition but evolved in 
two different frozen-N systems diverge from one another. 

2. Lyapunov exponents probe the rate of divergence asso- 
ciated with orbits evolved from the same initial condition 
in both a frozen-N system and in the smooth potential. 

This means that Lyapunov exponents also provide in- 
formation about the degree to which characteristics gen- 
erated in the smooth potential can be interpreted as pro- 
viding a pointwise approximation to real frozen- TV orbits 
with the same initial condition, as well as the degree to 
which orbits in two different frozen-TV systems remain 
close to one another in a pointwise sense. 

This is important at a practical level. In a real exper- 
iment one may perhaps be able to ensure that a given 
TV-body system constitutes (nearly) a fair sampling of 
some specified density distribution, but the details of the 
actual TV-body distribution are likely inaccessible. Of 



obvious interest, therefore, are the questions: to what 
extent will orbits in two different TV-body systems coin- 
cide? and to what extent do such orbits coincide with 
characteristics in the bulk potential associated with the 
smooth density distribution? 

Figure 16 exhibits the analogue of Figure 14, now 
generated by comparing the same 50 initial conditions 
evolved in two different frozen-TV systems. Figure 17 
compares orbits in a frozen- TV system with orbits in the 
smooth potential. In each case, the duration of the initial 
interval of especially fast exponential divergence is signif- 
icantly reduced, but the second interval with divergence 
at a rate ~ \S is still conspicuous. 

It is not hard to understand why the smooth poten- 
tial xs provides information about orbits in different 
frozen- TV simulations and/or their relation to orbits in 
the smooth potential. As noted already, discreteness ef- 
fects can be extremely well-mimicked by noise, at least 
mesoscopically. However, after the rapid decay of any 
initial transients, multiple noisy realisations of the same 
initial condition corresponding to a chaotic orbit typi- 
cally diverge exponentially in such a fashion that 24] 

SZ cx {Qn) 1 ' 2 eMxst) cx TV" 1 / 2 exp( X st), (36) 

where the second proportionality follows from the ob- 
served scaling rj cx 1/iV. 

V. CONCLUSIONS 

Viewed macroscopically, there is a precise sense in 
which, as TV increases, trajectories in frozen-TV sys- 
tems converge towards characteristics in the correspond- 
ing smooth potential. For very small particle number, 
TV < 10 4 or so, the notion of an average bulk potential 
fails and orbits in frozen-TV systems are very different 
from smooth potential characteristics. In particular, the 
usual distinctions between regularity and chaos that exist 
in a smooth potential seem completely lost. [28| However, 
for larger TV one begins to observe clear distinctions be- 
tween orbits evolved from initial conditions which, in the 
continuum limit, correspond to regular versus chaotic or- 
bits. 

In particular, although discreteness effects cannot be 
neglected, phase mixing of initially localised orbit ensem- 
bles in a frozen- TV environment allows for clear distinc- 
tions between 'regular' and 'chaotic' clumps. Just as for 
clumps evolved in a smooth potential, emittance growth 
for a regular frozen- TV clump proceeds as a power law 
in time, whereas it is roughly exponential for a chaotic 
clump. However, in both cases the growth is more rapid 
than in the smooth potential. Discreteness effects ac- 
celerate emittance growth for both regular and chaotic 
clumps. 

In terms of both the statistics of orbit ensembles and 
the complexity of individual orbits, discreteness effects 
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can be extremely well-mimicked by Gaussian white noise 
in the context of a Fokker-Planck, or Langevin, descrip- 
tion, with a coefficient of dynamical friction r\ and a 
diffusion constant D consistent with the predicted scal- 
ing D oc 77 oc (In A)/N, with In A the Coulomb logarithm. 
A Fokker-Planck/ Langevin description appears justified 
even when considering the short time behaviour of in- 
dividual orbits. This suggests strongly that Langevin 
simulations can be used to assess the importance of dis- 
creteness effects in systems where N is too large to allow 
honest direct summation integrations. 

To the extent that such an extrapolation is justified, 
one concludes that discreteness effects can remain im- 
portant even for very large N, especially for the case of 
chaotic orbits. Consider, e.g., the role of discreteness 
effects in accelerating emittance growth for an initially 
localised clump. For the case of the thermal equilibrium 
model, a relatively benign system without particularly 
large density contrasts and without internal substruc- 
tures, one needs N = 10 8 or more to justify the con- 
tinuum approximation in tracking the evolution of a reg- 
ular clump confined initially to a region ~ 1CP 3 the size 
of the accessible phase space. For the case of a chaotic 
clump of comparable size, one needs at least N = 10 9 5 . 
This has obvious implications for beams in that it affects 
macroscopic mixing and associated changes in the overall 
phase-space volume. Discreteness effects are also impor- 
tant because they can trigger transitions between regular 
and chaotic behaviour, a potentially serious problem for 
charged particle beams. One might, e.g.,, try to initialise 
a bunch in such a fashion that, although the bulk poten- 
tial admits chaotic orbits, only regular regions are popu- 
lated, thus aiming to facilitate emittance compensation. 
The problem, however, is that discreteness effects could 
transform significant numbers of orbits from regular to 
chaotic, thus making compensation far more difficult. 

The time scale associated with transitions between reg- 
ularity and chaos increases with increasing N, such tran- 
sitions being impossible in the continuum limit; but for 
any finite N there is presumably a maximum time over 
which it is safe to ignore these transitions. The criti- 
cal point, then, as is evident from Fig. 10, is that that 
time can be much shorter than the collisional relaxation 
time tu- To the extent that discreteness effects in the 
thermal equilibrium model can be mimicked by Gaus- 
sian white noise, particle number N — 10 6 corresponds 
to 77 ~ 10~ 4 ' 5 , which in turn implies a relaxation time 

~ V" 1 ~ 10 4 5 . It is, however, evident that, for a 
N = 10 6 realisation of the thermal equilibrium model, 
transitions from regularity to chaos can be important al- 
ready within a time t < 10 2 5 or so! By contrast, the 
dynals time scale tr> ~ 20. 

It should also be stressed that, even if discreteness ef- 
fects are too weak to facilitate frequent transitions be- 
tween regularity and chaos, they could well play an im- 
portant role in accelerating diffusion through a complex 



chaotic phase space. Generic smooth potentials admit- 
ting both regular and chaotic orbits have chaotic phase 
space regions partitioned by complex structures associ- 
ated with cantori in two dimensions and the Arnold web 
in three which, albeit not acting as absolute obstructions, 
serve as 'entropy' barriers to slow phase space trans- 
port |30| . The important point, then, is that even very 
low amplitude Gaussian white noise can dramatically ac- 
celerate diffusion through such barriers To the ex- 
tent that discreteness effects can be modeled as Gaussian 
white noise, they too should act as a significant source of 
accelerated phase space transport. 

The meaning of 'chaos' in the A~-body problem is nec- 
essarily somewhat subtle. In particular, it is important 
to recognise that two distinct sources of chaos can ex- 
ist, associated with physics on different scales. Short 
range microchaos, associated with close encounters be- 
tween individual charges, is a generic feature of the N- 
body problem, independent of the form of the bulk po- 
tential. However, there is also the possibility of larger 
scale macrochaos which arises if and only if, in the con- 
tinuum limit, the bulk potential admits chaos. The im- 
portant point, then, is that these two distinct sources of 
chaos can be characterised separately by different sets of 
Lyapunov exponents. Close encounters trigger an expo- 
nential separation of nearby trajectories at a rate xn- 
The bulk potential triggers an exponential separation at 
a rate xs which is typically much smaller. 

Standard numerical computations of Lyapunov expo- 
nents yield estimates of the much larger xn 7 a quantity 
that does not decrease with increasing N. This leads to 
the seemingly oxymoronic conclusion that the A^-body 
problem remains strongly chaotic for very large N, even 
if the potential is integrable in the N 00 limit and even 
if orbits 'look' nearly regular and have Fourier spectra 
that are nearly periodic. The crucial point, however, is 
that even though microchaos remains strong in the sense 
that xn does not decrease with increasing N, it becomes 
progressively less important macroscopically. The scale 
on which the exponential divergence saturates is compa- 
rable to a typical interparticle separation r sep , a distance 
that decreases as A^ -1 / 3 with increasing N. By track- 
ing the divergence of nearby orbits, starting with initial 
separations <C r sep and continuing until the separation 
becomes macroscopic, it is possible to extract estimates 
of both xn and xs- 

Finally, it should be noted that, as applied to the N- 
body problem, the smooth potential Lyapunov exponent 
Xs does not simply quantify the average divergence of 
two nearby trajectories in a single frozen-A^ simulation. 
It also quantifies the rate at which a frozen- N trajectory 
will diverge from a smooth potential characteristic with 
the same initial condition and the rate at which orbits 
with the same initial condition diverge in different frozen- 
N simulations, two quantities which, in some settings, 
could be even more important physically. 
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FIG. 1: The x and y coordinates of 1600 initially localised 
points evolved in frozen- N realisations of the regular potential 
(4) with N = 10 3 - 5 (left column) and N = W 5 (right col- 
umn) at different times t. From top to bottom, t = 16.0, 
t = 32.0, t = 64.0, t = 128.0, and t = 256.0. In each 
case, the integrations were performed with softening parame- 
ter e = 10~ 5 . 
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FIG. 2: The x and y coordinates of 1600 initially localised 
points evolved in frozen- iV realisations of the chaotic potential 
(5) with variable N at different times t. From left to right, 
one has N = 10 3 5 , N = 10 4 ' 5 , N = 10 5 5 and the smooth 
potential. From top to bottom, t — 16.0, t — 32.0, t = 64.0, 
t = 128.0, and t = 256.0. Once again e = 10~ 5 . 
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FIG. 3: The three-dimensional emittance 
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computed for the same regular clump used to generate FIG- 
URE 1, allowing for both frozen- iV backgrounds and energy- 
conserving white noise, (a) N = 10 3 . (b) 77 = 1(T 2 ' 5 (c) 



N = 10". (d) T) = 10" 3 (e) N = 10 4 . (f) r, 
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FIG. 4: The three-dimensional emittance e = {t x t y e z ) X ^' i 
computed for the same chaotic clump used to generate FIG- 
URE 2, allowing for both frozen-iV backgrounds and energy- 
conserving white noise. Note the logarithmic scale, (a) 



(c) N = 10 3 °. (d) 77 = 10" 



(e) 



N = 10 2 - 5 . (b) 77 = 10" 

N = 10 4 ' 5 . (f) 77 = 10" 40 . (g) N = 10 5 - 5 . (h) 77 = 10" 5 . (i) 
The smooth potential. The dashed line corresponds to a slope 
equaling the mean Lyapunov exponent (xs) for the orbits, (j) 
77 = 10~ 7 ' 5 , the strongest noise without an appreciable effect 
on e. 




FIG. 5: The x — y projection of a frozen-^ orbit generated 
from an initial condition corresponding in the smooth poten- 
tial to a circular orbit, (a) N = 10 2 5 . (b) N = 10 3 (c) 
N = 10 3 - 5 . (d) N = W 4 ° (e) N = 10 4 ' 5 . (f) N = 10 5 (g) 
N = 10 5 - 5 . (h) The smooth potential orbit. 
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FIG. 6: The power spectrum |x(w)| for the orbits exhibited 



in the preceding Figure, (a) N = 10 
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FIG. 7: (a) Diamonds show the complexity n(0.95), defined 
as the mean number of frequencies required to capture 95% of 
the total power, computed for a collection of 100 initial con- 
ditions integrated in frozen-iV realisations of the integrable 
Model 1 with variable N. Triangles show the same quantity 
for the same initial conditions integrated in the smooth po- 
tential but subjected to Gaussian white noise with variable rj 
The solid line exhibits the mean complexity for orbits evolved 
in the unperturbed smooth potential, (b) The same as (a), 
generated from the same initial conditions but now computed 
for the strongly chaotic Model 2. 
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FIG. 8: The three-dimensional emittance e computed for a 
clump of regular initial conditions for the thermal equilibrium 
model, allowing for both frozen-iV backgrounds and energy 
conserving white noise, (a) N = 10 4,J . (b) rj = 1CF 3 ' . (c) 
N = 10 5 - . (d) n = 1(T 3 ' 5 . (e) n = 10 5 ' 5 . (f) r) = 10" 40 . (g) 
N = 10 6 . (h) r) = 10 -4 - 5 . (i) Unperturbed motion in the 
smooth potential, (j) r\ = 1CF 6 ' 5 , the largest value of r\ that 
does not significantly impact emittance growth. 
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FIG. 9: The three-dimensional emittance e computed for a 
clump of chaotic initial conditions for the thermal equilibrium 
model, again allowing for both frozen- JV backgrounds and en- 
ergy conserving white noise, (a) N = 10 4 ' 5 . (b) r] = 10 _3 '°. 

10 5 - . (d) n = 1(T 3 ' 5 . (e) N = 10 5 ' 5 . (f) T) = 1(T 4 - . 
r,6 -° (h) t/ = 1CT 4 - 5 . (i) Unperturbed motion in the 
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smooth potential, (j) r) = 10 , the largest value of t] that 
does not significantly impact emittance growth. 
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FIG. 10: The percentage of frozen- iV orbits generated from a 
clump of regular initial conditions and evolved in the thermal 
equilibrium model which, at time t, have been converted to 
chaotic orbits. From top to bottom, the curves correspond 
to frozen- iV backgrounds with N = 10 4 ' 5 , N = 10 5 , and 
N = 10 5 ' 5 , and N = 10 6 () . 




FIG. 11: Mean value of the largest Lyapunov exponent \n 
as a function of iV for different choices of softening parameter 
e, computed for regular initial conditions evolved in Model 1. 
e = 1CF 5 : solid line and diamonds, e = 1CF 4 : dashed line and 
squares, e = 10~ 3 : dot-dashed line and triangles, e = 1CF 2 : 
triple-dot-dashed line and crosses, e = 10 _1 : dots and stars. 
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FIG. 12: Mean value of the largest Lyapunov exponent \n 
as a function of softening parameter e for different choices of 
N, computed for regular initial conditions evolved in Model 
1. N = 10 5 : solid line and pluses. N = 10 4 ' 5 : broad dashed 
line and stars. N — 10 4 : triple-dot-dashed line and crosses. 
N = 10 3 ' 5 : dot-dashed line and squares. N = 10 30 : dashed 
line and triangles. N = 10 2 ' 5 : dots and diamonds. 
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FIG. 13: Mean value of the largest Lyapunov exponent \n 
as a function of N for the same initial condition integrated in 
the integrable Model 1 (solid curve) and the chaotic Model 2 
(short dashed curve). In both cases, e = 10~°. The broad- 
dashed line at the bottom corresponds to the smooth poten- 
tial Lyapunov exponent \s for the same initial condition in- 
tegrated in the smooth potential corresponding to Model 2. 
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FIG. 14: The mean phase space separation |<5Z(i)| for 50 
nearby pairs of initial conditions evolved in frozen-iV reali- 
sations of the chaotic Model 2 with (from top to bottom) 
N = 10 5 , 10 5 ' 5 , and 10 6 . The solid line has a slope 
Xs = 0.056, equal to the the mean value of the largest smooth 
potential Lyapunov exponent. The dashed line has a slope 
Xjv = 0.82, equal to the mean value of the largest JV-body 
exponent. 




FIG. 15: The same as the preceding for regular orbits in the 
thermal equilibrium model. The dashed line again has a slope 
equal to the mean value of the largest iV-body Lyapunov ex- 
ponent. Note the absence of the second exponential phase. 
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FIG. 16: The mean phase space separation |<5Z(i)| for 50 
initial conditions evolved in two different frozen- N realisations 
of Model 2 with (from top to bottom) N = 10 5 , 10 55 , and 
10 6 . The solid line again has a slope \s = 0.056, equal 
to the mean value of the largest smooth potential Lyapunov 
exponent. 




FIG. 17: The mean phase space separation |(5Z(i)| for 50 
initial conditions evolved in Model 2, both in the smooth po- 
tential and in a frozen- N system with (from top to bottom) 
N = 10 5 , 10 5 ' 5 , and 10 6 . The solid line again has a slope 
Xs = 0.056, equal to the mean value of the largest smooth 
potential Lyapunov exponent. 



